Diversity and the nature of responses to environmental change.
Owen Petchey & Jennifer Probst
9/28/2018
1 Setup R
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(deSolve)
library(rootSolve)
set.seed(13)
run_sims <- FALSE2 Introduction
The aim is to reproduce the model based results / figures.
Here is the original publication: https://www.nature.com/articles/s41467-017-00912-x
And here is the supplementary information, including the table of parameter values: https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-017-00912-x/MediaObjects/41467_2017_912_MOESM1_ESM.pdf
Note that below I go with the notation from the table S1 and not the main text. However there is only one small difference between these: In the main text the notation is simplified to y_S_PB (which is y_SR_PB in Table S1) and y_S_SB (which is y_SO_SB in Table S1).
The published model is modified:
- Random noise in substrate concentrations can be included.
- Temporal changes in oxygen diffusitivity can be included.
- Immigration can be included.
- A minimum abundance of each species can be set.
3 Methods
3.1 Model dimensions
- Time: hours
- Volume: litres
- Substrate quantity: micromoles
- Organism quantity: cells
3.2 Dynamic model
source("functions/bushplus_dynamic_model.r")3.3 Defaults
3.3.1 Default parameter values
default_parameter_values <- c(
## maximum specific growth rates
g_max_CB = 0.05,
g_max_PB = 0.07,
g_max_SB = 0.1,
## half-saturation constants
k_PB_SR = 10,
k_SB_SO = 5,
k_CB_P = 0.2,
k_PB_P = 0.5,
k_SB_P = 0.5,
## half-inhibition constant
h_SR_CB = 300,
h_O_PB = 100,
h_O_SB = 100,
## yield (cells per micromole)
Y_SO_SB = 3.33e7,
Y_SR_PB = 1.25e7, ## Y_S_PB in the main text
y_P_CB = 1.67e8,
y_P_PB = 1.67e8,
y_P_SB = 1.67e8,
## oxygen production per microbial cell
P_CB = 6e-9,
## mortality rate
m_CB = 0.02,
m_PB = 0.028,
m_SB = 0.04,
## substrate diffusivity
a_S = 0.001,
a_O = 8e-4,
a_P = 0.01,
## background substrate concentration
back_SR = 300,
back_SO = 300,
back_O = 300,
back_P = 9.5,
## immigration rates (cells per time)
i_CB <- 0,
i_PB <- 0,
i_SB <- 0,
## oxidisation rate of reduced sulphur
c = 4e-5
)3.3.2 Default event
3.3.2.1 Default event definition
Used to add noise and to ensure minimum abundance. Noise is added to the resource concentrations, and not to the organism abundances. Minimum organism abundance is set. I.e. organismal abundance cannot go below some very low level.
source("functions/default_event_definition.r")3.3.2.2 Default event schedule
This is when an event happens
default_event_interval <- 103.3.2.3 Default noise event options
Default to no noise:
default_noise_sigma <- 03.3.2.4 Default minimum abundance event
Default to following minimums
default_minimum_abundances <- rep(1e2, 3)
names(default_minimum_abundances) <- c("CB", "PB", "SB") 3.3.3 Default simulation duration
default_sim_duration <- 2000
default_sim_sample_interval <- 13.3.4 Default oxygen diffusivity (a) forcing
Here we specify no forcing of a_0 (used unless set elsewhere).
default_log10a_series <- c(log10(default_parameter_values["a_O"]),
log10(default_parameter_values["a_O"]))3.3.5 Default initial state
Initial conditions favouring anoxic
default_initial_state <- c(N_CB = 5e1,
N_PB = 1e7,
N_SB = 1e7,
SO = 300,
SR = 300,
O = 1e1,
P = 1e1)3.4 Simulation function
source("functions/run_simulation.r")3.5 Helper functions
source("functions/plot_dynamics.r")
source("functions/ss_expt_a_N.r")
source("functions/get_nonlinearity.r")
source("functions/get_hysteresis_range.r")3.6 Code test
# dynamic_model = bushplus_dynamic_model
# event_definition = default_event_definition
# parameter_values = default_parameter_values
# event_interval = default_event_interval
# noise_sigma = default_noise_sigma
# minimum_abundances = default_minimum_abundances
# sim_duration = default_sim_duration
# sim_sample_interval = default_sim_sample_interval
# log10a_forcing = default_log10a_forcing
# initial_state = default_initial_state
# solver_method = "radau"
#
#
# times <- seq(0,
# sim_duration,
# by = sim_sample_interval)
#
# event_times <- seq(min(times),
# max(times),
# by = event_interval)
#
# log10a_forcing <- approxfun(x = log10a_forcing[,1],
# y = log10a_forcing[,2],
# method = "linear",
# rule = 2)
# assign("log10a_forcing", log10a_forcing, envir = .GlobalEnv)
#
# out <- as.data.frame(ode(y = initial_state,
# times = times,
# func = dynamic_model,
# parms = parameter_values,
# method = solver_method,
# events = list(func=event_definition,
# time=event_times)))
#
# sim_res <- run_simulation(dynamic_model = bushplus_dynamic_model,
# event_definition = default_event_definition,
# parameter_values = default_parameter_values,
# event_interval = default_event_interval,
# noise_sigma = default_noise_sigma,
# minimum_abundances = default_minimum_abundances,
# sim_duration = default_sim_duration,
# sim_sample_interval = default_sim_sample_interval,
# log10a_forcing = default_log10a_forcing,
# initial_state = default_initial_state,
# solver_method = "radau")
#sim_res <- run_simulation()
#plot_dynamics(sim_res)4 Results: Reproducing previous results
4.1 Initial conditions favouring anoxic, no noise, no forcing
Figure 2a, b in Bush et al.
init_state_favouring_anoxic <- c(N_CB = 5e1,
N_PB = 1e7,
N_SB = 1e7,
SO = 300,
SR = 300,
O = 1e1,
P = 1e1)
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic)
plot_dynamics(simulation_result = sim_res)
#ggsave('bacteria_no_noise.png')Population dynamics seem as in the paper. Resource dynamics fluctuate more than in the paper.
4.2 Initial conditions favouring anoxic, with noise, no forcing
Figure 2a, b in Bush et al, but with noise
## add some noise
noise_sigma_0.01 <- 0.01
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
noise_sigma = noise_sigma_0.01)
plot_dynamics(sim_res)
#ggsave("high_SB_PB.png")4.3 Initial conditions favouring oxic, no noise, no forcing
Figure 2c, d in Bush et al.
init_state_favouring_oxic <- c(N_CB = 1e7,
N_PB = 1e5,
N_SB = 1e5,
SO = 500,
SR = 50,
O = 30,
P = 10)
sim_res <- run_simulation(initial_state = init_state_favouring_oxic)
plot_dynamics(sim_res)
##ggsave("high_CB.png")Population and resource dynamics seen as those in the paper.
4.4 Initial conditions favouring oxic, with noise, no forcing
Figure 2c, d in Bush et al. but with noise.
sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
noise_sigma = noise_sigma_0.01)
plot_dynamics(sim_res)
##ggsave("high_CB.png")5 Results: Composition/diversity and nature of response
ssfind_minimum_abundances <- rep(0, 3)
names(ssfind_minimum_abundances) <- c("CB", "PB", "SB")
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
# ssfind_init_state <- c(N_CB = NA,
# N_PB = NA,
# N_SB = NA,
# SO = 500,
# SR = 50,
# O = 30,
# P = 10)
##ssfind_init_state <- default_initial_state
ssfind_init_state <- init_state_favouring_anoxic
ssfind_init_state <- c(N_CB = NA, #this one is varied
N_PB = 1e8,
N_SB = 1e8,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
grid_num_a <- 100
a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
grid_num_N <- 25.0.1 Each composition
5.0.1.1 No organisms
initial_CBs <- 0
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_no_biology <- ss_by_a_N(expt)
saveRDS(ss_no_biology, "data/ss_no_biology.RDS")ss_no_biology <- readRDS("data/ss_no_biology.RDS")ss_no_biology %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_no_biology %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.2 CB
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_CB <- ss_by_a_N(expt)
saveRDS(ss_only_CB, "data/ss_only_CB.RDS")ss_only_CB <- readRDS("data/ss_only_CB.RDS")tempx1 <- ss_only_CB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
geom_path(size=0.8) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
theme_light() +
ylim(0,10) +
geom_path(data = filter(tempx1, Species == "PB"),
aes(y=0.2), size = 0.8) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="Only CB") 
ggsave("figures/CB.png", width = 3, height = 2)ss_only_CB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.3 SB
initial_CBs <- 0 # this one is varied
initial_PBs <- 0
initial_SBs <- 10^seq(0, 8, length=grid_num_N)
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_SB <- ss_by_a_N(expt)
saveRDS(ss_only_SB, "data/ss_only_SB.RDS")ss_only_SB <- readRDS("data/ss_only_SB.RDS")ss_only_SB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.4 PB
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_only_PB <- ss_by_a_N(expt)
saveRDS(ss_only_PB, "data/ss_only_PB.RDS")ss_only_PB <- readRDS("data/ss_only_PB.RDS")tempx1 <- ss_only_PB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.))
tempx1 %>%
ggplot(aes(x = a, y = log10(Quantity), col = Species)) +
geom_path(size = 0.8) +
theme_light() +
ylim(0,10) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
geom_path(data = filter(tempx1, Species == "CB"),
aes(y=0.2), size = 0.8)+
geom_path(data = filter(tempx1, Species == "PB"),
size = 0.8) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="Only PB") 
ggsave("figures/PB.png", width = 3, height = 2)ss_only_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.5 SB and PB
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_SB <- expt$N_PB
ss_only_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_PB, "data/ss_only_SB_PB.RDS")ss_only_SB_PB <- readRDS("data/ss_only_SB_PB.RDS")ss_only_SB_PB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.6 SB and CB
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 0
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_SB <- expt$N_CB
ss_only_SB_CB <- ss_by_a_N(expt)
saveRDS(ss_only_SB_CB, "data/ss_only_SB_CB.RDS")ss_only_SB_CB <- readRDS("data/ss_only_SB_CB.RDS")ss_only_SB_CB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_CB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:5) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.7 CB and PB
initial_CBs <- 0 # this one is varied
initial_PBs <- 10^seq(0, 8, length=grid_num_N)
initial_SBs <- 0
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
expt$N_CB <- expt$N_PB
ss_only_CB_PB <- ss_by_a_N(expt)
saveRDS(ss_only_CB_PB, "data/ss_only_CB_PB.RDS")ss_only_CB_PB <- readRDS("data/ss_only_CB_PB.RDS")ss_only_CB_PB %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_only_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.1.8 CB, SB, and PB
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
ss_CB_SB_PB <- ss_by_a_N(expt)
saveRDS(ss_CB_SB_PB, "data/ss_CB_SB_PB.RDS")ss_CB_SB_PB <- readRDS("data/ss_CB_SB_PB.RDS")ss_CB_SB_PB %>%
select(a, CB=N_CB, SB=N_SB, PB=N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species),
size=0.8) +
ylab(expression(Log[10](Abundance))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
theme_light() +
ylim(0,10) +
theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="CB, SB, & PB") 
ggsave("figures/CB_SB_PB.png", width = 3, height = 2)ss_CB_SB_PB %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
5.0.2 All response measures
all_ss <- bind_cols(readRDS("data/ss_no_biology.RDS"),
composition = rep("No organisms", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a)))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_CB.RDS"),
composition = rep("CB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB.RDS"),
composition = rep("SB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_PB.RDS"),
composition = rep("PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB_PB.RDS"),
composition = rep("SB-PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_SB_CB.RDS"),
composition = rep("SB-CB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_only_CB_PB.RDS"),
composition = rep("CB-PB", grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
all_ss <- bind_rows(all_ss,
bind_cols(readRDS("data/ss_CB_SB_PB.RDS"),
composition = rep("CB-SB-PB",
grid_num_a*grid_num_N),
direction = c(rep("up", grid_num_a),
rep("down", grid_num_a))))
x_long <- all_ss %>%
gather(key = Species, value = Quantity, 2:8)
x_wide <- x_long %>%
select(-initial_N_CB) %>%
spread(key = direction, value=Quantity, drop=T)
resp_summary <- x_wide %>%
group_by(composition, Species) %>%
summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
hyst_2 = get_hysteresis_range(up, down, a),
nl_up = get_nonlinearity(a, log10(up)),
nl_down = get_nonlinearity(a, log10(down))) %>%
ungroup() %>%
gather(key = resp_measure, value = value, 3:6) %>%
mutate(composition = fct_relevel(composition,
"No organisms",
"CB",
"SB",
"PB",
"CB-PB",
"SB-CB",
"SB-PB",
"CB-SB-PB"),
num_func_groups = case_when(composition =="No organisms" ~ 0,
composition == "CB" ~ 1,
composition == "SB" ~ 1,
composition == "PB" ~ 1,
composition == "CB-PB" ~ 2,
composition == "SB-CB" ~ 2,
composition == "SB-PB" ~ 2,
composition == "CB-SB-PB" ~ 3
),
var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
)
saveRDS(resp_summary, "data/resp_summary.RDS")resp_summary <- readRDS("data/resp_summary.RDS")
resp_summary %>%
ggplot() +
geom_point(aes(x = composition, y = value)) +
facet_grid(Species ~ resp_measure, scales = "free") +
coord_flip()
5.0.3 Nonlinearity
5.0.3.1 All species nonlinearity
resp_summ1 <- resp_summary %>%
filter(resp_measure %in% c("nl_down", "nl_up")) %>%
group_by(composition, Species, num_func_groups, var_type) %>%
summarise(ave_nl = mean(value))
resp_summ1 %>%
ggplot() +
geom_point(aes(composition, ave_nl)) +
facet_wrap( ~ Species)
5.0.3.2 Substrate nonlinearity
focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(composition) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
5.0.3.3 Organism nonlinearity
focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(composition) %>%
summarise(ave_ave_nl = sum(ave_nl)/unique(num_func_groups))
ave_nl_organism %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))## Warning: Removed 1 rows containing missing values (geom_point).
5.0.3.4 System nonlinearity
focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
group_by(composition) %>%
summarise(ave_ave_nl = sum(ave_nl) / (unique(num_func_groups)+4))
ave_nl_system %>%
ggplot() +
geom_point(aes(composition, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
ggsave("figures/system_nonlinearity.png", width = 3, height = 2)5.0.4 Hysteresis
5.0.4.1 Substrate hysteresis
focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(composition, num_func_groups, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = composition, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
5.0.4.2 Organism hysteresis
focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(composition, num_func_groups, resp_measure, var_type) %>%
summarise(ave_hyst = sum(value) / unique(num_func_groups))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = composition, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))## Warning: Removed 2 rows containing missing values (geom_point).
5.0.4.3 System hysteresis (separate)
resp_summ4 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
group_by(resp_measure, composition) %>%
summarise(mean_hyst = sum(value) / (unique(num_func_groups)+4))
resp_summ4 %>%
ggplot() +
geom_point(aes(x = composition, y = mean_hyst)) +
facet_wrap( ~ resp_measure) +
ggtitle("Hysteresis in system response.")
5.0.4.4 System hysteresis (combined)
resp_summ5 <- resp_summ4 %>%
group_by(composition) %>%
summarise(mean_mean_hyst = mean(mean_hyst))
resp_summ5 %>%
ggplot() +
geom_point(aes(x = composition, y = mean_mean_hyst)) +
xlab("Functional groups present") +
ylab("Amount of hysteresis") +
ylim(0,2) +
theme_light() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("figures/composition_system_hysteresis.png", width = 2.5, height = 2.5)6 Results: Varying strength of inhibition
6.1 Introduction
Three parameters set the strength of inhibition (half-inhibition constants):
- h_SR_CB = 300,
- h_O_PB = 100,
- h_O_SB = 100,
These work by the inhibition function:
\(\frac{1}{1 + \frac{X}{h_X}}\)
inhibition## function (X, h_X)
## 1/(1 + (X/h_X))Plot some examples:
default_parameter_values["h_SR_CB"]## h_SR_CB
## 300SR <- 10^seq(1.4, 2.5, length = 100)
ggplot() +
geom_line(aes(x = SR, y = inhibition(SR, default_parameter_values["h_SR_CB"]))) +
ylab("Multiplicative effect\non CB growth rate")
default_parameter_values["h_O_PB"]## h_O_PB
## 100O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_PB"]))) +
ylab("Multiplicative effect\non PB growth rate")
default_parameter_values["h_O_SB"]## h_O_SB
## 100O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, default_parameter_values["h_O_SB"]))) +
ylab("Multiplicative effect\non SB growth rate")
proportion_change <- 2
default_parameter_values["h_O_SB"]## h_O_SB
## 100O <- 10^seq(-0.2, 2.5, length = 100)
ggplot() +
geom_line(aes(x = O, y = inhibition(O, proportion_change*default_parameter_values["h_O_SB"]))) +
ylab("Multiplicative effect\non SB growth rate")
6.2 Experiment
Increasing the inhibition half saturation constants weakens the inhibition.
ssfind_minimum_abundances <- default_minimum_abundances
ssfind_simulation_duration <- 10000
ssfind_simulation_sampling_interval <- ssfind_simulation_duration
ssfind_event_interval <- ssfind_simulation_duration
ssfind_parameters <- default_parameter_values
ssfind_init_state <- c(N_CB = NA, #this one is varied
N_PB = NA,
N_SB = NA,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
grid_num_mult <- 20
grid_num_N <- 2
grid_num_a <- 100
a_Os <- 10^seq(-4.5, -1.5, length=grid_num_a)
inhib_const_mults <- 2^seq(-1, 0.4, length=grid_num_mult)
initial_CBs <- 10^seq(0, 8, length=grid_num_N)
initial_PBs <- 1e5
initial_SBs <- 1e5
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
inhib_const_mult = inhib_const_mults,
a_O = a_Os)ss_h1 <- ss_by_a_h_N(expt)
saveRDS(ss_h1, "data/ss_h1.RDS")ss_h1 <- readRDS("data/ss_h1.RDS")6.3 Results
6.3.1 Close to default inhibition
inhib_const_mults## [1] 0.5000000 0.5262004 0.5537737 0.5827919 0.6133306 0.6454696 0.6792927
## [8] 0.7148882 0.7523489 0.7917725 0.8332620 0.8769256 0.9228772 0.9712366
## [15] 1.0221302 1.0756906 1.1320576 1.1913783 1.2538074 1.3195079focal_inhibition_multiplier <- inhib_const_mults[14]
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
6.3.2 Double strength inhibition
## YES! 0.5 does mean twice as strong inhibition
focal_inhibition_multiplier <- 0.5
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
6.3.3 Weaker inhibition
## YES! 2 does mean half as strong inhibition
focal_inhibition_multiplier <- 2^0.4
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, N_CB, N_SB, N_PB) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
ss_h1 %>%
filter(inhib_const_mult == focal_inhibition_multiplier) %>%
select(a, SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 2:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species))
6.3.4 Nonlinearity
x_long <- ss_h1 %>%
mutate(direction = ifelse(initial_N_CB==1, "up", "down")) %>%
gather(key = Species, value = Quantity, 3:9)
x_wide <- x_long %>%
select(-initial_N_CB) %>%
spread(key = direction, value=Quantity, drop=T)
resp_summary <- x_wide %>%
group_by(inhib_const_mult, Species) %>%
summarise(hyst_1 = mean(abs(log10(up) - log10(down))),
hyst_2 = get_hysteresis_range(up, down, a),
nl_up = get_nonlinearity(a, log10(up)),
nl_down = get_nonlinearity(a, log10(down))) %>%
ungroup() %>%
gather(key = resp_measure, value = value, 3:6) %>%
mutate(
var_type = ifelse(str_sub(Species, 1, 1)=="N", "Organism", "Substrate")
)6.3.4.1 All species nonlinearity
resp_summ1 <- resp_summary %>%
filter(resp_measure %in% c("nl_down", "nl_up")) %>%
group_by(inhib_const_mult, Species, var_type) %>%
summarise(ave_nl = mean(value))
resp_summ1 %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_nl)) +
facet_wrap( ~ Species)
6.3.4.2 Substrate nonlinearity
focus_on <- "Substrate"
ave_nl_substrate <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_substrate %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
6.3.4.3 Organism nonlinearity
focus_on <- "Organism"
ave_nl_organism <- resp_summ1 %>%
filter(var_type == focus_on) %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_organism %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
6.3.4.4 System nonlinearity
focus_on <- "System"
ave_nl_system <- resp_summ1 %>%
group_by(inhib_const_mult) %>%
summarise(ave_ave_nl = mean(ave_nl))
ave_nl_system %>%
ggplot() +
geom_point(aes(inhib_const_mult, ave_ave_nl)) +
ggtitle(paste("Nonlinearity in", focus_on, "response."))
ggsave("figures/system_nonlinearity.png", width = 3, height = 2)6.3.5 Hysteresis
6.3.5.1 Substrate hysteresis
focus_on <- "Substrate"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(inhib_const_mult, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
6.3.5.2 Organism hysteresis
focus_on <- "Organism"
resp_summ3 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2"),
var_type == focus_on) %>%
group_by(inhib_const_mult, resp_measure, var_type) %>%
summarise(ave_hyst = mean(value))
resp_summ3 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = ave_hyst)) +
facet_wrap( ~ resp_measure, scales="free") +
ggtitle(paste("Hysteresis in", focus_on, "response."))
6.3.5.3 System hysteresis (separate)
resp_summ4 <- resp_summary %>%
filter(resp_measure %in% c("hyst_1", "hyst_2")) %>%
group_by(resp_measure, inhib_const_mult) %>%
summarise(mean_hyst = mean(value))
resp_summ4 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = mean_hyst)) +
facet_wrap( ~ resp_measure) +
ggtitle("Hysteresis in system response.")
6.3.5.4 System hysteresis (combined)
resp_summ5 <- resp_summ4 %>%
group_by(inhib_const_mult) %>%
summarise(mean_mean_hyst = mean(mean_hyst))
resp_summ5 %>%
ggplot() +
geom_point(aes(x = inhib_const_mult, y = mean_mean_hyst)) +
xlab("Organismal tolerance") +
ylab("Amount of hysteresis") +
theme_light() +
ylim(0,2)
ggsave("figures/inhibition_system_hysteresis.png", width = 2.5, height = 2)7 Results: Temporal state switching
7.1 Switch anoxic-oxic-anoxic
## anoxic initial state
init_state <- init_state_favouring_anoxic
## simulation timing
simulation_duration <- 5000
simulation_sampling_interval <- 10
event_interval_t <- 100
## low-high-low oxgygen diffusivity
lhl_log10a_series <- c(-8, -8, -2, -2, -2, -2, -2, -8, -8)
## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_anoxic,
log10a_series = lhl_log10a_series,
sim_duration = simulation_duration,
sim_sample_interval = simulation_sampling_interval,
event_interval = event_interval_t)
saveRDS(sim_res, "data/lhl_with_defaults.RDS")lhl_with_defaults <- readRDS("data/lhl_with_defaults.RDS")
plot_dynamics(lhl_with_defaults)
lhl_with_defaults$result %>%
filter(between(time, 1000, 73000)) %>%
ggplot() +
geom_path(aes(x=a, y=log10(O)),
arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))7.2 Switch oxic-anoxic-oxic
## anoxic initial state
init_state <- init_state_favouring_oxic
## simulation timing
simulation_duration <- 50000
simulation_sampling_interval <- 10
event_interval_t <- 100
## low-high-low oxgygen diffusivity
hlh_log10a_series <- c(-2, -2, -10, -14, -18, -22, -2, -2, -2)
## run simulation
sim_res <- run_simulation(initial_state = init_state_favouring_oxic,
log10a_series = hlh_log10a_series,
sim_duration = simulation_duration,
sim_sample_interval = simulation_sampling_interval,
event_interval = event_interval_t)
saveRDS(sim_res, "data/hlh_with_defaults.RDS")hlh_with_defaults <- readRDS("data/hlh_with_defaults.RDS")
plot_dynamics(hlh_with_defaults)
hlh_with_defaults$result %>%
filter(between(time, 1000, 100000)) %>%
ggplot() +
geom_path(aes(x=a, y=log10(O)),
arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")))8 Appendices
8.1 runsteady function not finding correct values
init_state_anoxic_2 <- c(N_CB = 1e1,
N_PB = 1e8,
N_SB = 1e8,
SO = 300,
SR = 300,
O = 1e-5,
P = 1e1)
RS <- runsteady(y = init_state_anoxic_2,
fun = bushplus_dynamic_model,
parms = default_parameter_values,
times = c(0, 1e5),
stol=1)
RS$y
RS <- runsteady(y = init_state_favouring_oxic,
fun = bushplus_dynamic_model,
parms = default_parameter_values,
times = c(0, 1e5),
stol=1)
RS$y8.2 Separatrix
Note that this section not working / updated since the code was made into functions
8.2.1 Effect of N_CB and oxygen diffusability
noise_sigma <- 0.0
get_final_states_N_CB <- function(x) {
#give N_CB, N_SB, N_PB and a_O values and put these into state and parameters
state["N_CB"] <- x["N_CB"]
parameters["a_O"] <- x["a_O"]
## update forcing function
a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
log10(parameters["a_O"]),
max(times),
log10(parameters["a_O"])))
log10a_forcing2 <- approxfun(x = a_forcing1[,1],
y = a_forcing1[,2],
method = "linear",
rule = 2)
## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)
as.data.frame(ode(y = state,
times = times,
func = micro_mod_var,
parms = parameters,
method = "radau",
events = list(func=noise_events, time=times)))[2,-1]
#gives you values for the last time step
}
times <- seq(0, 50000, by = 50000)
state <- c(N_CB = NA, #this one is varied
N_PB = 1e8,
N_SB = 1e8,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
#different combinations of N_CB and oxygen diffusability
grid_num <- 50
expt <- expand.grid(N_CB = 10^seq(0, 8, length=grid_num),
a_O = 10^seq(-3.5, -2, length=grid_num))
result <- apply(expt, 1, function(x) get_final_states_N_CB(x)) %>%
tibble() %>%
unnest() %>%
mutate(initial_N_CB=expt$N_CB,
a_O=expt$a_O,
highlow=ifelse(N_CB>2e8, 'high', 'low'))
result <- result %>%
mutate(N_CB = ifelse(N_CB < minimum_abundance, minimum_abundance, N_CB))
result$highlow <- as.factor(result$highlow)
saveRDS(result, "data/CB_oxdiff_hyst.RDS")result <- readRDS("data/CB_oxdiff_hyst.RDS")#plot N_CB / P_CB against a_0
p <- ggplot(result) +
geom_line(aes(x=log10(a_O), y=log10(N_CB), color=highlow), size=1.25) +
ggtitle("Stable states of the system")+
xlab('log[oxygen diffusivity]') +
ylab('Initial cyanobacteria abundance')
p <- p + geom_point(data=result, aes(x=log10(a_O), y=log10(initial_N_CB), color=highlow)) + labs(color = "Final CB concentration")
plot(p)
#ggsave("initial conditions_1.png")#show hysterisis effects for all the parameters
#plot oxygen against diffusability
ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(O))))+
xlab('log[oxygen diffusivity]') +
ylab('log[oxygen concentration]')+
geom_segment(aes(x=log10(result$a_O[839]), y=log10(result$O[839]), xend=log10(result$a_O[838]), yend=log10(result$O[838])), arrow=arrow(), colour="blue", size=0.5)+
geom_segment(aes(x=log10(result$a_O[1761]), y=log10(result$O[1761]), xend=log10(result$a_O[1762]), yend=log10(result$O[1762])), arrow=arrow(), colour="blue", size=0.5)+
geom_segment(aes(x=-4.25, y=-0.45, xend=-3.75, yend=0.1), arrow=arrow(), colour="blue", size=0.5)+
geom_curve(aes(x = -3.75, y = 1.9, xend = -4.25, yend = 1.5), arrow=arrow(), colour='blue', size = 0.5, curvature = -0.2)
#ggsave("hysterisis_CB_1.png")ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(P))))+
ggtitle("Hysterisis effect of phoshate concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[phosphate concentration]')
#ggsave("hysterisis_CB_2.png")ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(SO))))+
ggtitle("Hysterisis effect of oxidized sulfur concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[oxydized sulfur concentration]')
#ggsave("hysterisis_CB_3.png")ggplot(result) +
geom_point(aes(x=log10(a_O), y=(log10(SR))))+
ggtitle("Hysterisis effect of reduced sulfur concentration")+
xlab('log[oxygen diffusivity]') +
ylab('log[reduced sulfur] concentration]')
#ggsave("hysterisis_CB_4.png")These graphs show an at least qualitative match to the result in the paper.
8.2.2 Effect of PB/SB depending on oxygen diffusability
noise_sigma <- 0.00
get_final_states_N_PB <- function(x) {
#give in N_CB, N_SB, N_PB and a_O values and put these into state and parameters
state["N_PB"] <- x["var"]
state["N_SB"] <- x["var"]
parameters["a_O"] <- x["a_O"]
## update forcing function
a_forcing1 <- matrix(ncol=2, byrow=T, data=c(0,
log10(parameters["a_O"]),
max(times),
log10(parameters["a_O"])))
log10a_forcing2 <- approxfun(x = a_forcing1[,1],
y = a_forcing1[,2],
method = "linear",
rule = 2)
## assign the changed a_forcing2 into the global environment, from where the model gets a_forcing2
assign("log10a_forcing2", log10a_forcing2, envir = .GlobalEnv)
as.data.frame(ode(y = state,
times = times,
func = micro_mod_var,
parms = parameters,
method = "radau", events = list(func=noise_events, time=times)))[2,-1]
#gives you values for the last time step
}
times <- seq(0, 50000, by = 5000)
state <- c(N_CB = 5,
N_PB = NA,#this one is varied
N_SB = NA,
SO = 200,
SR = 200,
O = 10,
P = 9.5)
#different combinations of N_PB and oxygen diffusivity
grid_num <- 50
expt <- expand.grid(var=10^seq(0, 7, length=grid_num),
a_O = 10^seq(-3, -2, length=grid_num))
result2 <- apply(expt, 1, function(x) get_final_states_N_PB(x)) %>%
tibble() %>%
unnest() %>%
mutate(initial_N_PB=expt$var,
initial_N_SB=expt$var,
a_O=expt$a_O,
highlow=ifelse(N_PB>1e8, 'high', 'low'))
result2 <- result2 %>%
mutate(N_PB = ifelse(N_PB < minimum_abundance, minimum_abundance, N_PB))
result2$highlow <- as.factor(result2$highlow)
saveRDS(result2, "data/PBSB_oxdiff_hyst.RDS")result2 <- readRDS("data/PBSB_oxdiff_hyst.RDS")#plot N_PB against a_0
p <- ggplot(result2) +
geom_line(aes(x=log10(a_O), y=log10(N_PB), color=highlow), size=1.25) +
ggtitle("Stable states of the system")+
xlab('log[oxygen diffusivity]') +
ylab('Initial PS and SR bacteria concentration')
p <- p + geom_point(data=result2,
aes(x=log10(a_O), y=log10(initial_N_PB), color=highlow))
p <- p + labs(color = "Final PB and SB concentration")
plot(p)
#ggsave("initial conditions_2.png")8.3 EWS analyses
Evaluation of the switch: Get the functions:
source('C:/Users/probs/Desktop/Sem 6/Internship/functions_analysis_EWS_2.R')Get state shifts and calculate parameters:
x3 <- 30 #size of wished rolling window (percentage btw 0 and 100)
x5 <- 50 #size of wished rolling window (percentage btw 0 and 100)
#first switch
name_1='anoxic_oxic 1'
trans1 <- data.frame(out_t %>% filter(time>87500 & time<92500)) #for analysis
trans1_f <- data.frame(out_f %>% filter(time>87500 & time<92500)) #for plotting
n_13 <- round(length(trans1$time)*x3/100) #size of rolling window
t_13 <- length(trans1$time) - n_13 + 1 #number of rolling windows
n_15 <- round(length(trans1$time)*x5/100) #size of rolling window
t_15 <- length(trans1$time) - n_15 + 1 #number of rolling windows
out_df_1_N_CB <- data.frame(trans1$N_CB)
colnames(out_df_1_N_CB) <- c("N_CB")
out_df_1_N_PB <- data.frame(trans1$N_PB)
colnames(out_df_1_N_PB) <- c("N_PB")
out_df_1_N_SB <- data.frame(trans1$N_SB)
colnames(out_df_1_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans1_f) +
ggtitle("Anoxic_oxic 1")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 1.png")
#second switch
name_2='oxic_anoxic 1'
trans2 <- data.frame(out_t %>% filter(time>335000 & time<360000))
trans2_f <- data.frame(out_f %>% filter(time>335000 & time<360000))
n_23 <- round(length(trans2$time)*x3/100)
t_23 <- length(trans2$time) - n_23 + 1
n_25 <- round(length(trans2$time)*x5/100)
t_25 <- length(trans2$time) - n_25 + 1
out_df_2_N_CB <- data.frame(trans2$N_CB)
colnames(out_df_2_N_CB) <- c("N_CB")
out_df_2_N_PB <- data.frame(trans2$N_PB)
colnames(out_df_2_N_PB) <- c("N_PB")
out_df_2_N_SB <- data.frame(trans2$N_SB)
colnames(out_df_2_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans2_f) +
ggtitle("Oxic_anoxic 1")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 1.png")
#third switch
name_3='oxic_anoxic 2'
trans3 <- data.frame(out2_t %>% filter(time>165000 & time<185000))
trans3_f <- data.frame(out2_f %>% filter(time>165000 & time<185000))
n_33 <- round(length(trans3$time)*x3/100)
t_33 <- length(trans3$time) - n_33 + 1
n_35 <- round(length(trans3$time)*x5/100)
t_35 <- length(trans3$time) - n_35 + 1
out2_df_1_N_CB <- data.frame(trans3$N_CB)
colnames(out2_df_1_N_CB) <- c("N_CB")
out2_df_1_N_SB <- data.frame(trans3$N_SB)
colnames(out2_df_1_N_SB) <- c("N_SB")
out2_df_1_N_PB <- data.frame(trans3$N_PB)
colnames(out2_df_1_N_PB) <- c("N_PB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans3_f) +
ggtitle("Oxic_anoxic 2")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("oxic_anoxic 2.png")
#fourth switch
name_4='anoxic_oxic 2'
trans4 <- data.frame(out2_t %>% filter(time>340000 & time<348500))
trans4_f <- data.frame(out2_f %>% filter(time>340000 & time<348500))
n_43 <- round(length(trans4$time)*x3/100)
t_43 <- length(trans4$time) - n_43 + 1
n_45 <- round(length(trans4$time)*x5/100)
t_45 <- length(trans4$time) - n_45 + 1
out2_df_2_N_CB <- data.frame(trans4$N_CB)
colnames(out2_df_2_N_CB) <- c("N_CB")
out2_df_2_N_PB <- data.frame(trans4$N_PB)
colnames(out2_df_2_N_PB) <- c("N_PB")
out2_df_2_N_SB <- data.frame(trans4$N_SB)
colnames(out2_df_2_N_SB) <- c("N_SB")
ggplot(aes(x=time, y=trans_quantity, col=species), data=trans4_f) +
ggtitle("Anoxic_oxic 2")+
xlab('Time') +
ylab('Trans_quantity') +
facet_wrap(~var_type, scales = "free", ncol = 1) + #create one plot for organisms and one for substrate
geom_line()+
labs(color = "Organism/Substrate")
ggsave("anoxic_oxic 2.png")Rolling window analysis
#analysis for 1st
rollingwindow_analysis(trans1, 2, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 3, n_15, t_15, name_1)
rollingwindow_analysis(trans1, 4, n_15, t_15, name_1)
#analysis for 2nd
rollingwindow_analysis(trans2, 2, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 3, n_25, t_25, name_2)
rollingwindow_analysis(trans2, 4, n_25, t_25, name_2)
#analysis for 3rd
rollingwindow_analysis(trans3, 2, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 3, n_35, t_35, name_3)
rollingwindow_analysis(trans3, 4, n_35, t_35, name_3)
#analysis for 4th
rollingwindow_analysis(trans4, 2, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 3, n_45, t_45, name_4)
rollingwindow_analysis(trans4, 4, n_45, t_45, name_4)Metric based indicators from the earlywarnings package
#analysis for 1st
earlywarnings_plot(out_df_1_N_CB, name_1, x3)
earlywarnings_plot(out_df_1_N_PB, name_1, x3)
earlywarnings_plot(out_df_1_N_SB, name_1, x3)
#analysis for 2nd
earlywarnings_plot(out_df_2_N_CB, name_2, x3)
earlywarnings_plot(out_df_2_N_PB, name_2, x3)
earlywarnings_plot(out_df_2_N_SB, name_2, x3)
#analysis for 3rd
earlywarnings_plot(out2_df_1_N_CB, name_3, x3)
earlywarnings_plot(out2_df_1_N_PB, name_3, x3)
earlywarnings_plot(out2_df_1_N_SB, name_3, x3)
#analysis for 4th
earlywarnings_plot(out2_df_2_N_CB, name_4, x3)
earlywarnings_plot(out2_df_2_N_PB, name_4, x3)
earlywarnings_plot(out2_df_2_N_SB, name_4, x3)Conditional Heteroskedasticity
#analysis for 1st
ch <- condhet_analysis(out_df_1_N_CB, name_1)
ch <- condhet_analysis(out_df_1_N_PB, name_1)
ch <- condhet_analysis(out_df_1_N_SB, name_1)
#analysis for 2nd
ch2 <- condhet_analysis(out_df_2_N_CB, name_2)
ch2 <- condhet_analysis(out_df_2_N_PB, name_2)
ch2 <- condhet_analysis(out_df_2_N_SB, name_2)
#analysis for 3rd
ch3 <- condhet_analysis(out2_df_1_N_CB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_PB, name_3)
ch3 <- condhet_analysis(out2_df_1_N_SB, name_3)
#analysis for 4th
ch4 <- condhet_analysis(out2_df_2_N_CB,name_4)
ch4 <- condhet_analysis(out2_df_2_N_PB, name_4)
ch4 <- condhet_analysis(out2_df_2_N_SB, name_4)Detrending fluctuation analysis
#analysis for 1st
dfa_analysis(trans1, 2, n_15, t_15, name_1)
dfa_analysis(trans1, 3, n_15, t_15, name_1)
dfa_analysis(trans1, 4, n_15, t_15, name_1)
#analysis for 2nd
dfa_analysis(trans2, 2, n_25, t_25, name_2)
dfa_analysis(trans2, 3, n_25, t_25, name_2)
dfa_analysis(trans2, 4, n_25, t_25, name_2)
#analysis for 3rd
dfa_analysis(trans3, 2, n_35, t_35, name_3)
dfa_analysis(trans3, 3, n_35, t_35, name_3)
dfa_analysis(trans3, 4, n_35, t_35, name_3)
#analysis for 4th
dfa_analysis(trans4, 2, n_45, t_45, name_4)
dfa_analysis(trans4, 3, n_45, t_45, name_4)
dfa_analysis(trans4, 4, n_45, t_45, name_4)Drift Diffusion Jump Nonparametrics Early Warning Signals:
#shift 1
ddj11<- ddj_model(out_df_1_N_CB)
plot_ddj_1(out_df_1_N_CB, ddj11, name_1)
plot_ddj_2(out_df_1_N_CB, ddj11, name_1)
ddj12<- ddj_model(out_df_1_N_PB)
plot_ddj_1(out_df_1_N_PB, ddj12, name_1)
plot_ddj_2(out_df_1_N_PB, ddj12, name_1)
ddj13<- ddj_model(out_df_1_N_SB)
plot_ddj_1(out_df_1_N_SB, ddj13, name_1)
plot_ddj_2(out_df_1_N_SB, ddj13, name_1)
#shift 2
ddj21<- ddj_model(out_df_2_N_CB)
plot_ddj_1(out_df_2_N_CB, ddj21, name_2)
plot_ddj_2(out_df_2_N_CB, ddj21, name_2)
ddj22<- ddj_model(out_df_2_N_PB)
plot_ddj_1(out_df_2_N_PB, ddj22, name_2)
plot_ddj_2(out_df_2_N_PB, ddj22, name_2)
ddj23<- ddj_model(out_df_2_N_SB)
plot_ddj_1(out_df_2_N_SB, ddj23, name_2)
plot_ddj_2(out_df_2_N_SB, ddj23, name_2)
#shift 3
ddj31<- ddj_model(out2_df_1_N_CB)
plot_ddj_1(out2_df_1_N_CB, ddj31, name_3)
plot_ddj_2(out2_df_1_N_CB, ddj31, name_3)
ddj32<- ddj_model(out2_df_1_N_PB)
plot_ddj_1(out2_df_1_N_PB, ddj32, name_3)
plot_ddj_2(out2_df_1_N_PB, ddj32, name_3)
ddj33<- ddj_model(out2_df_1_N_SB)
plot_ddj_1(out2_df_1_N_SB, ddj33, name_3)
plot_ddj_2(out2_df_1_N_SB, ddj33, name_3)
#shift 4
ddj41<- ddj_model(out2_df_2_N_CB)
plot_ddj_1(out2_df_2_N_CB, ddj41, name_4)
plot_ddj_2(out2_df_2_N_CB, ddj41, name_4)
ddj42<- ddj_model(out2_df_2_N_PB)
plot_ddj_1(out2_df_2_N_PB, ddj42, name_4)
plot_ddj_2(out2_df_2_N_PB, ddj42, name_4)
ddj43<- ddj_model(out2_df_2_N_SB)
plot_ddj_1(out2_df_2_N_SB, ddj43, name_4)
plot_ddj_2(out2_df_2_N_SB, ddj43, name_4)